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Estimating the spectral characteristics of a nonstationary random process is an important but challenging 

Cn , task, which can be facilitated by exploiting structural properties of the process. In certain applications, the 

observed processes are underspread, i.e., their time and frequency correlations exhibit a reasonably fast decay, 

and approximately time-frequency sparse, i.e., a reasonably large percentage of the spectral values is small. 

For this class of processes, we propose a compressive estimator of the discrete Rihaczek spectrum (RS). 

d \ This estimator combines a minimum variance unbiased estimator of the RS (which is a smoothed Rihaczek 

C/^^ \ distribution using an appropriately designed smoothing kernel) with a compressed sensing technique that exploits 

the approximate time-frequency sparsity. As a result of the compression stage, the number of measurements 

^ ^ required for good estimation performance can be significantly reduced. The measurements are values of the 

*^ ' discrete ambiguity function of the observed signal at randomly chosen time and frequency lag positions. We 

tH- ' provide bounds on the mean-square estimation error of both the minimum variance unbiased RS estimator and 

the compressive RS estimator, and we demonstrate the performance of the compressive estimator by means of 



in 



^ZJ ' simulation results. The proposed compressive RS estimator can also be used for estimating other time-dependent 

P\J ■ spectra (e.g., the Wigner-Ville spectrum) since for an underspread process most spectra are almost equal. 

^^ Index Terms 

"X 

\^ , Nonstationary spectral estimation, time-dependent power spectrum, Rihaczek spectrum, Wigner-Ville spee- 



ch 



trum, compressed sensing, basis pursuit, cognitive radio. 

I. Introduction 

Estimating the spectral characteristics of a random process is an important task in many signal analysis 
and processing problems. Conventional spectral estimation based on the power spectral density is restricted 
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to wide-sense stationary and, by extension, wide-sense cyclostationary processes lH], ||2]. However, in many 
applications — including speech and audio, communications, image processing, computer vision, biomedical 
engineering, and machine monitoring — the signals of interest cannot be well modeled as wide-sense (cy- 
clo)stationary processes. In particular, in cognitive radio systems Q-llSl, the receiver has to infer from the 
received signal the location of unoccupied frequency bands ("spectral holes") that can be used for data 
transmission. Here, modeling the received signal as a nonstationary process can be advantageous because 
it potentially allows a faster estimation of time- varying changes in band occupation 131. 

For a general nonstationary process, a "power spectral density" that is nonnegative and extends all the 
essential properties of the conventional power spectral density is not available |[6l- |[T0l . Several different 
definitions of a "time-dependent (or time-varying) power spectrum" have been proposed in the literature, 
see |[6l- |[26l and references therein. A time-dependent power spectrum is a function of both frequency and 
time. Examples include the Wigner-Ville spectrum ||2]-||9], US, Ell, HTl, Rihaczek spectrum (RS) Q, IS, 
El, (21, evolutionary spectrum ||T3], ||22], (23], ||29], Weyl spectinm (23l, Page specti-um E], (H), Levin 
spectrum (T2l . and physical spectrum lITSl . ll30l . Some of these spectra are not guaranteed to be nonnegative 
or even real-valued; in general, nonnegativity is paid for by the loss of certain other desirable properties. 
However, it has been shown (TOl . |[24l that in the practically important case of nonstationary processes with 
fast decaying time-frequency (TF) correlations — so-called underspread processes (TOl . (24ll - (26l . (3T]-(34] — 
all major spectra yield effectively identical results, are (at least approximately) real-valued and nonnegative, 
and satisfy several other desirable properties at least approximately. Thus, in the underspread case, the specific 
choice of a spectrum is of secondary theoretical importance and can hence be guided by practical considerations 
such as computational complexity. 

Once a specific definition of time-dependent spectrum has been adopted, an important problem is the 
estimation of the spectrum from a single observed realization of the process. This nonstationary spectral 
estimation problem is fundamentally more difficult than spectral estimation in the (cyclo)stationary case, 
because long-term averaging cannot be used to reduce the mean-square error of the estimate. Formally, any 
estimator of a nonparametric time-dependent spectrum can also be viewed as a TF representation of the 
observed signal (7], (261 . (351 . (361 . Estimators have been previously proposed for several spectra including 
the Wigner-Ville specti-um and RS (e.g., (2]-(9], ED, (26], El, ITTl-Bll). 

In this paper, extending our work in ll42l . we propose a "compressive" estimator of the RS that uses the 
recently inttoduced methodology of compressed sensing (CS) (43l . ll44l . The proposed estimator is suited 
to underspread processes that are approximately TF sparse. The latter property means that only a moderate 
percentage of the values of the discrete RS are significantly nonzero. Both assumptions — underspreadness and 



TF sparsity — are reasonably well satisfied in typical cognitive radio applications. We consider the RS because it 
is the simplest time-dependent spectrum from a computational viewpoint, especially in the discrete setting used. 
The proposed compressive estimator of the RS is obtained by augmenting a basic noncompressive estimator 
(a smoothed version of the Rihaczek distribution, cf. Q, Ol, 1261, 1361, Ell, 61, ED, gl-ETl) with a 
CS compression-reconstruction stage. 

Compressive spectral estimation methods have been proposed previously, also in the context of cognitive 
radio Il48l . ||49l . However, these methods are restricted to the estimation of the power spectral density of 
stationary processes. Also, they perform CS directly with respect to the observed signal (process realization), 
whereas our method performs CS with respect to an estimate of a TF autocorrelation function known as the 
expected ambiguity function (EAF). This EAF estimate is an intermediate step in the calculation of the spectral 
estimator, somewhat similar to a sufficient statistic. In some sense, we perform a twofold compression, first 
by using only an EAF estimate (instead of the raw observed data) for spectral estimation and secondly by 
"compressing" that estimate. This approach can be advantageous if dedicated hardware units for computing 
values of the EAF estimate are employed, because fewer such units are required. It can also be advantageous 
if the values of the EAF estimate have to be transmitted over low-rate links — e.g., in wireless sensor networks 
ll50l — or stored in a memory, because fewer such values need to be transmitted or stored. 

A major focus of our work is an analysis of the estimation accuracy of the proposed compressive estimator. 
Because finding a closed-form expression of the mean-square error is intractable, we derive upper bounds on 
the mean-square error. These bounds depend on two components: the first component is determined by the 
amount of "underspreadness," corresponding to the concentration of the EAF of the observed process; the 
second component is related to the TF sparsity properties of the observed process. As we will see below, there 
is a tradeoff between these components, since a well concentrated EAF of an underspread process tends to 
imply a poorly concentrated RS, which is disadvantageous in terms of TF sparsity. 

The remainder of this paper is organized as follows. In Section |lll we state our general setting and review 
some fundamentals of nonstationary random processes and their TF representation. In Section |III] we describe 
a basic noncompressive estimator of the RS. In Section |IVl we develop a compressive estimator by augmenting 
the noncompressive estimator with a CS compression-reconstruction stage. Bounds on the mean-square error 
of both the noncompressive and compressive estimators are derived in Section [V] Finally, numerical results 
are presented in Section |Vl] 

Notation. The modulus, complex conjugate, real part, and imaginary part of a complex number a G C are 
denoted by \a\, a*, R{a}, and ^{a}, respectively. Boldface lowercase letters denote column vectors belonging 
to C^ for some L G N, whereas boldface uppercase letters denote matrices belonging to C*^^^ for some 



M, N £ N. The kth entry of a vector a is denoted by (a)^, and the entry of a matrix A in the ith row and jth 
column by (A)^ . The superscripts ^ and ^ denote the transpose and Hermitian transpose, respectively, of a 
vector or matrix. The ii norm of a vector a G C^ is denoted by ||a||^ = X]fc=i I(^)a:I' ^^'^ ^^^ ^2 norm by 



||a||2 = VsL^a. The number of nonzero entries is denoted by ||a||Q. The trace of a square matrix A G ([^MxM 
is denoted by tr{A} = Xlfcii {■^)k,k- Given a matrix A G C*^^^, we denote by vec{A} G C*^^ the vector 
obtained by stacking all columns of A. Given two matrices A G C*^^^^^ and B G C*^^^^^ we denote by 
A(g)B G C^'^i^'i^x^i^^ their Kronecker product [51]. The inner product of two square matrices A, B G C*^^^^ 
is defined as (A, B) = tr{AB^}. The Kronecker delta is denoted by 6[m\, i.e., 6[m\ G {0, 1} and ^[tti] = 1 if 
and only if m = 0. Finally, [N] = {0,1,. ..,N- 1}. 

II. Expected Ambiguity Function and Rihaczek Spectrum 

In this section, we state our setting and review some fundamentals of the TF representation of nonstationary 
random processes. Let Xjn], n£ [N\ be a nonstationary, zero-mean, circularly symmetric complex, finite-length 
or, equivalently, periodic^ random process with autocorrelation function (ACF) ^x[ni,n2] = E{X[ni]X*[n2]}, 
where E{-} denotes expectation. Since we observe process samples only for nG [N], we consider the ACF 
7x[ni,n2] only for ni,n2 G [N]. This is justified for a process that is well concentrated in the interval 
[N]. An equivalent representation of the ACF ^x[ni,n2] is the correlation matrix Vx — E{xx^}, where 
X ^ (X[0] X[l] ■ ■ ■ X[iV-l])^G C^; note that (rx)„^+i,„^+i = lx[ni,n2] for ni,n2 G [N]. 

We assume that X[n] is an underspread process lITOl . Il24l - |[26l . |[3TI - |[34l . which means that its correlation 
in time and frequency decays reasonably fast. The underspread property is phrased mathematically in terms 
of the discrete EAF, which is defined as the following discrete Fourier transform (DFT) of the ACF lITOl . 

isi-na, EB, El, m- 

Ax[m,l] ^ Yl lx[n,n-m]^e-^'^'^. (1) 

ne[Af] 

Here, m and I denote discrete time lag and discrete frequency lag, respectively, and [ni, 722]^ = [ni mod N, 722 
mod A^]. Note that this definition of Ax[m, I] is iV -periodic in both m and I. The EAF Ax[m, I] is a TF-lag 
representation of the second-order statistics of X[n\ that describes the TF correlation structure of X[n]. A 
nonstationary process X[n\ is said to be underspread if its EAF is well concentrated around the origin in the 
(m, /)-plane, i.e., 

Ax[m,l]^0, \/{m,l)^A, with A={-M,...,M}j^x{-L,...,L}j^, 



where < M < 



N 



, 0<L< 



N 



, and ML < N. (2) 



It will be convenient to consider length- A'^ functions as periodic functions with period A'^. 



Here, e.g., {— M, . . . , M}j^ denotes the A^ -periodic continuation of the interval {— M, . . . , M}. The concen- 
tration of the EAF around the origin can be measured by the EAF moment defined in Section IV-AI (see (|43])). 
For later reference, we note that the EAF is the expectation of the ambiguity function (AF) 171, 1351 . |[36l 

Ax[m,l]^ J]X[n]X*[n-m]^e-^'^'", (3) 

n(i[N] 

i.e., Ax[m,l] =E{^xK/]}. 

Nonstationary spectral estimation is the problem of estimating a "time-dependent power spectrum" of 
the nonstationary process X[n] from a single realization x[n] observed for n G [N]. As mentioned earlier, 
there is no definition of a "time-dependent power spectrum" that satisfies all desirable properties ll6l- llT0l . 
However, in the underspread case considered, most reasonable definitions of a time-dependent power spectrum 
are approximately equal, represent the mean energy distribution of the process over time and frequency, and 
approximately satisfy all desirable properties lITOl . lE4l . We therefore use the simplest such definition, which 
is the RS Q, H, mi, m. The discrete RS is defined as the following DFT of the ACE: 

i?x[n,A;] ^ J]7x[n,n-m]^e-^^i^'='". (4) 

Just as the EAF and AF, the RS Rxln, k] is A^-periodic in both its variables. Furthermore, the RS is complex- 
valued in general, but it is approximately real-valued and nonnegative in the underspread case lITOl . Il24l . 
Averaging the RS over the discrete frequency variable k yields the mean power of X[n\, i.e., 

1 Y, Rx[nM = 7x[n,n] = EJI^Hp}. 

k(i[N] 

The RS is related to the EAF via a symplectic two-dimensional (2D) DFT: 

Rx[nM = ^ E ^xKZle-^'te^— "0, (5) 

m,l<^[N\ 

Ax[m,l\ = ^ J] i?x[n,fe]e^'i^(™'^-'"). (6) 

The relation dSjl extends the Fourier transform relation between the power spectral density and the autocorre- 
lation function of a stationary process lH], ||2] to the nonstationary case. It follows from (JS) that the RS of an 
underspread process is a smooth function. Furthermore, the RS is the expectation of the Rihaczek distribution 

(RD) m, m, m, Ea, eh 

Rx[n,k] ^ Y. X[n]^X*[n-m]^e-^'^^'^ = X[n]^X*[k]^e-^'^'^\ 



where X[k] = E„e[iv] ^N^"^^''" is the DFT of X[n]. That is, Rx[n,k\ = E{i?x[n, /c]}. The 2D DFT 
relation (jSj, (l6]l holds also for the RD and AF, i.e., 



(V) 



Ax[m,l] = ^ Y. ^x[n,fc]e^'^('-^-^"). 

n,ke[N] 

Our central assumption, besides the underspread property, is that the nonstationary process X[n] is "ap- 
proximately TF sparse" in the sense that only a moderate percentage of the RS values Rx[n,k] within the 
fundamental {n,k)-ve.gion [N^ = [N] x [A^] is significantly nonzero. For such approximately TF sparse 
processes, we will develop a compressive estimator of the RS by augmenting a basic RS estimator with a 
compression-reconstruction stage. We present the basic estimator first. 

III. Basic RS Estimator 

In analogy to well-known estimators of the Wigner-Ville spectrum C]-||9], EH, ll26l . BOl . 11411 . a basic 
(noncompressive) estimator of the RS Rx[n,k] is given by the following smoothed version of the RD ||28]| . 

ED: 

Rx[n,k] = ^ Y. ^n-n',k-k']Rx[n\k']. (8) 

n',k'e[N] 

Here, ^[n,k\ is a smoothing function that is A^-periodic in both arguments and whose choice will be discussed 
presently. We note that Rx[n, k] is a quadratic functional of the observed signal X[n], and it commutes with 
cyclic time and frequency shifts of the type X[n] ^^ e^~N^^X[n — rri\^ Q, |[36l . Because of Q, the symplectic 
2D inverse DFT of Rx[n, k], 

Ax[^,l\ =^ E i?x[n,A;]e^-t('"'=-'"), (9) 

n,ki^[N\ 

can be viewed as an estimator of the EAF Ax[m, I]. Using ([8]l and ([7]) in Q, we obtain 

Ax[m,l] = </>[m,/]Ax[m,/], (10) 

where the 2D window (weighting, taper) function (/>[m,/] is related to the smoothing function <l>[n, A;] through 
a 2D DFT, i.e., 

(f)[m,l] ^ — Y ^[n,A:]e^'i^("^-'"). (11) 

n,k&[N] 

Note that </>[m,/] and Ax[w,, /] are A^-periodic in both m and /. 



We now consider the choice of the smoothing function ^[n,k] or, equivalently, of the window function 
(/)[7n, /]. Our performance criterion is the mean-square error (MSE) 

e ^ E{\\Rx - RxWl) = Yl mRx[n,k]-Rx[n,kf}. 

n,k&[N] 

The MSE can be decomposed as s = B^ + V with the squared bias term B^ = \\E{Rx} — RxWo ^^id 
the variance V = E|||i?x — Eji^xllL}- We will consider a minimum variance unbiased (MVU) design of 
<I>[n,A;]. This means that Rx[n,k] is required to be unbiased, i.e., i? = 0, and the variance V is minimized 
under this constraint. More specifically, we will adopt the MVU design proposed in ||26l, ||4H, which is based 
on the idealizing assumption that the EAF Ax [m, I] is supported on a cyclically extended rectangular region 
A = {-M, ...,M}^ X {-L,...,L}jv, i.e., Ax[m,l] = for all {m,l) A, with < Af < [N/2\ and 
< L < \_N/2\. This is somewhat similar to the underspread property (0; however, it is an exact, rather than 
approximate, support constraint. As a further difference from the underspread property, we do not require that 
ML <^ N. We note that this idealizing exact support constraint is only needed for the MVU interpretation of 
our design of <l>[n, k]; in particular, it will not be used for our performance analysis in Section |Vl The size of 
A — i.e., the choice of L and M — is a design parameter that can be chosen freely in principle. The resulting 
estimator Rx,M"^v[n, k] (cf. (IT6l )) can be applied to any process X[n], including, in particular, processes whose 
EAF ^x[^n^;^] is not exactly supported on A. 

We briefly review the derivation of the MVU smoothing function presented in ll26l . ll4n . Using (ITOl ) and 
E{Ax['m,l]] = Ax[m,l], the bias term B'^ = \\E{Rx} — Rx\\2 = ||E{^x} — ^x||2 can be expressed as 

B^ = Y. Mm,l]-l)Ax[m,lf. (12) 

m,l&[N] 

Thus, i?^ = if and only if (j)[m, /] = 1 on the support of Ax [m, I], i.e., for all {m, I) G A. Under the constraint 
B^ = 0, minimizing the variance of Rx [n, k] is equivalent to minimizing the mean power 

P ^ E{\\Rx\\l} i Ejllixll^} ^E{\\cl,[m,l]Ax[m,l]\\l} = ^ \^[m,l]\^E{\Ax[m,lf}. 

Splitting this sum into a sum over [N] D A (where </)[m, I] = 1) and a sum over [N] D A (here, A denotes 
the complement of A), it is clear that P is minimized if and only if the latter sum is zero, which means that 
(l)[m, I] must be zero for {m, I) G [N] DA, and further, due to the periodicity of (l)[m, I], for (m, /) G A. Thus, 
we conclude that the MVU window function (DFT of the MVU smoothing function) is the indicator function 



/4[m, /] of the EAF support A = {-M, ..., M}^ x {-L, ... , L}^: 



ct>MWv[m,l\ = lA[m,l\ = < (13) 



A 

0, otherwise. 



The corresponding EAF estimator in (|9]l, (ITO^ is obtained as 

{Ax\rn I] (m l)£jA 
(14) 
0, otherwise. 

Therefore, the MVU estimator of the RS is given by (see ©) 



Rx,MMn,k] = ^ E ^x,Mvu[m,Z]e-^'t('=™-'^') 



(15) 



m,l&[N] 
M L 



m=-M l=-L 

where the periodicity of the summand with respect to m and I has been exploited in the last step. 

IV. Compressive RS Estimator 

Next, we will augment the basic RS estimator presented in the previous section with a compression- 
reconstruction stage. 

A. A Result from CS Theory 

The proposed compressive RS estimator is based on the following result from CS theory Il44l . Il53l . Let 
U be a unitary (up to a factor c > 0) and equimodular matrix of size Q x Q, i.e.. 

Furthermore let the "measurement matrix" M of size P x Q, with P <Q, be formed by randomly selecting 

P rows from U. (In practice, typically, P<^Q.) Consider a vector r of length Q that is observed through M, 

i.e., the observed (measured) vector is 

m = Mr, (17) 

which has length P. If P < Q, the above "measurement equation" Mr = m is underdetermined, so r cannot 
be reconstructed unambiguously from m unless we use further constraints. Therefore, we take the sparsest 
solution, which is formulated, in an approximate sense, as the r' with minimum £i norm that satisfies the 
measurement equation: 



f = arginin ||r'||-|^ 

Mr' = m 



This is a linear program that is known as Basis Pursuit (BP). The following has been shown 1441 . |[53l : for 
any i^ G N, if the number of measurements is chosen as 

P>C{logQ)^K, (18) 



then the £2 error ||f — r||2 satisfies with overwhelming probabilit 



I 



D „ 

|r-r||2 < -=\\r-rK\\i. (19) 



Here, vk is the K-sparse approximation of r, which is obtained by zeroing all entries of r except the K 
entries with largest magnitudes; furthermore, C and D are two positive constants that do not depend on r. 
The message of ( fT9l ) is that r will be close to the true vector r if r is approximately K-sparse. 

We will need the following extension of ( fT9l ). Let ^ C {1, . . . , Q} be an arbitrary index set of size |^| = i^, 
and let r^ denote the vector that is obtained from r by zeroing all entries of r except the K entries whose 
indices are in ^. It can then be easily verified that 

||r-r;^||i < ||r-r^||,. (20) 

Combining (l20l ) and ( fT9l ) yields the following result: If P > C (log Q)^K, then with overwhelming probability 



D 



for any index set ^ C {1, . . . , Q} of size \Q\ =K. 

B. Basic DFT Relation 

The proposed compressive RS estimator is furthermore based on a 2-D DFT relation that will now be 
derived. We recall from (IT4] | that the EAF estimate Ax,Mvu["i5 ^] is exactly zero outside the effective EAF 
support A = {-M, ..., M}^ x {-L, ... , L}^, where < M < [N/2\ and < L < [N/2\ . In what follows, 

we will denote by 

S = |[iV]^n^| = (2M + 1)(2L+1) (22) 

the size of one period of A. Because 2M + 1 and 2L + 1 do not necessarily divide A^, we furthermore 
define an "extended effective EAF support" as the periodized rectangular region A' = {—M, . . . , —M + 
AM — l}jy X {— L, . . . , —L + AL — l}jy, where AM and AL are chosen as the smallest integers such that 
^That is, the probability of ( I19t not being true decreases exponentially with P. 



10 



AAf > 2M + 1 and AL > 2L + 1 and, moreover, AM and AL divide N, i.e, there are integers An, AA; 
such that An AL = A/c AM = A^ or, equivalently. 



The size of one period of A' is 



N N 

^™ = XT > ^^ = XT7 • (23) 

AL ' AM 



S' 4 |[iV]2n^'| = AMAL. 



Note that 

^C^', 5<5', (24) 

although typically S ^ S'. Let us arrange the values of one period of Ax^uwi'm', I] that are located within A' 
into a matrix A G C^^xAL^ j ^^ 

(A)^+i^,+i ^ ix,Mvu[m -M,l-L], m G [AM] , / G [AL] . (25) 

Alternatively, we can represent Ax,Myu['m,l] by the matrix R G C^^^^*^ whose entries are given by the 
following 2D DFT of dimension AL x AM: 

«^ > «^ > -o { q(m-M) p(l-L)\ 

me [AM] «e[AL] 
-Af+AA/-1 -L+AL-1 



*^ E E ix,Mvu[m,/]e-^-2'^(i^-^) 



m=-M l=-L 

M L 



Y^ YAx[m,l]e-^^^^^-^) , pG [AL], gG [AM]. (27) 



{HI 

m=-Ml=-L 



It can be seen by comparing (l27l ) and (fT6l ) that the matrix entries (R) _|_]^ ^^ equal (up to a constant factor) 
a subsampled version of Rx,Mvvi[n,k], i.e., 

Wp+i.g+i = A^iix,Mvu[pAn, gA/c] , p G [AL] , q G [AM] , (28) 

with An = N/ AL and AA; = N/AM as in (l23l) . This subsampUng does not cause a loss of information 
because Ax,mvu["i, /] is supported in A, and therefore, by (|24] |. also in ^' = {— M, . . . , — M + AM — 1}^ x 
{_L,...,_L + AL-1}^. 
Inverting (l26l ). we obtain 

(A)„+M+i = ^7 E E i^Ui,,+i^''''^'"^''^^^ mG [AM], ZG [AL]. (29) 

pe[AL]g6[AA/] 



11 

This 2-D DFT relation will constitute an important basis for our compressive RS estimator. It can be compactly 

written as 

U^r = a , (30) 

where r = vec{R} G C"^', a = vec{A} G C"^', and 

U = -^Ffi^FAAf e C^'^^'. (31) 

Here, e.g., Fal denotes the DFT matrix for DFT length AL, which is defined as 

(FALW,+i^e-^-^-"^, p,l€[AL]. 
Furthermore, using (l25i in (ITST l. we obtain 

me[AAf]Ze[AL] 

Inserting (l29l ). we see that the basic RS estimate iix,Mvu[^5 k] can be calculated from r (or, equivalently, from 
R) according to 

-Rx,MVu[?^,^] = C{r}[n,k] 






me [AM] is [AL] 



{m-M-)q (l-L)p 



pe[AL]q6[AAf] 



X e-J'i^[''(™-*^)-"('-^)] . (32) 



C. Measurement Equation and Sparse Reconstruction 

The compressive RS estimator can be obtained by combining the results of the previous two subsections. To 
motivate our development, we assume that i?x,Mvu[pAn, q/S.k] is approximately K-sparse for some K < S', i.e., 
at most K of the 5" values of the basic RS estimator i?x,Mvu[™> k] on the subsampled grid (n, k) = (pAn, qAk) 
are significantly nonzero. (Because Rx,M\i][n, k] is an estimator of the RS, this assumption is consistent with 
our basic assumption that the RS Rxinjk] itself is approximately sparse.) Due to (|28] ). it follows that the 
matrix R and, equivalently, the vector r = vec{R} are approximately K-sparse. Furthermore, according to 
(l30l) . r G C-^' is related to the EAF estimate a = vec{A} G C"^' as U^r = a, where U^ (see ^T^) is an 
orthogonal and equimodular matrix. Motivated by our discussion in Section ITV-A[ let us define a*-^) G C^ as 
the vector made up of P randomly selected entries of a, for some P < S' (typically, P<^ S'). Thus, recalling 
(|25] ) and (IT4l) . the entries of a*-^) are P values of the masked AF Iji^[m,l]Ax[m,l] randomly located within 



12 

the region [N] D A' or, equivalentlylj the values of /^[rri, IjAxim, I] at P randomly chosen TF lag positions 
(m, /) G {-M, ..., -M + AM - 1} x {-L, . . . , -L + AL - 1}. We then have from ^ 

Mr = a(^) , 

where the Px S' matrix M is obtained by randomly selecting P rows from U^; the indices of these rows 
correspond to the indices of the entries selected from a. 

The above equation is recognized to be an instance of the measurement equation (ITtI i. with m in dTTl ) 
given by a^^^. From our discussion in Section HV-AI we can then conclude the following: For 

P > C (log S')^K = C [ log(AM) + log(AL)] V (33) 

(see (dD), the result of BP operating on a^^^ i.e., 



f = argmin ||r||j^, (34) 

Mr' = a(-P) 



satisfies with high probability 



D 



.G\ 



^„ „i, (35) 



for any index set Q of size |^| =if (see (|2TI)). Here, as before, r^ denotes the vector that is obtained from r 
by zeroing all entries of r except the K entries whose indices are in Q. Since r is approximately A'-sparse, 
the index set Q can be chosen such that the corresponding entries {{^)k\keQ comprise, with high probabilityo 
the significantly nonzero entries of r, implying a small norm ||r — r^||^. The bound (1351 ) then shows that BP 
is capable of reconstructing r — and, thus, the subsampled basic RS estimator i?x,Mvu[pAn, (/A/;;] — from the 
compressed AF vector a^^^ with a small error. (We recall, at this point, that the entries of r equal the values 
of ^x,Mvub^"'i gAA;].) 

D. The Compressive RS Estimator 

From the BP reconstruction result f in (l34l ). a compressive approximation of the basic RS estimator 
Rx,MNv[^: /c] in (IT6l ) is finally obtained by substituting f for r in (l32l ): 

^Typically, the region [A'^] n A' is only slightly larger than the effective EAF support [TV] n .4. Thus, most of the P entries of a'^' 

are values of Ax [m, I] randomly located within [A*'] n A or, equivalently, {— M, . . . , M} x {—L, . . . , L}. The remaining entries of 

(P) 
a^ ' are zero. 

■*Note that the index set Q is deterministic and fixed, whereas the indices of the largest entries of r may vary with each realization. 
However, for the performance analysis in Section |V] it is sufficient to assume that the index set Q approximately contains the indices 
of the largest entries of r for each realization. 
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Rx,cs[n,k] = C{r}[n,k] 



]^ E E E E *u,„«»'=-( 

me[AM]Ze[AL] 



AM AL / 



pe[/\L] qe[AM] 

X e-Ji^['=(™-*^)-"('-^)] , (36) 

where R = unvecjr) is the matrix corresponding to f . This defines the compressive RD estimator. 

To summarize, the proposed compressive RD estimator Rx,cs[''T',k] is calculated by the following steps. 

1) Choose K such that it reflects the prior intuition about the effective sparsity of the subsampled RS 
Rx[p^n,qAk]. (Equivalently, KN'^/S' reflects the prior intuition about the effective sparsity of the RS 

Rx[n,k].) 

2) Acquire P > C[log(AM) + log(AL)] K values of the masked AF Ij([m,l] Ax[m,l] at randomly 
chosen TF lag positions (m, /) G {-M, . . . , -M + AM - 1} x {-L, . . . , -L + AL - 1}. Let a(^) 
denote the vector containing these "compressive measurements." A compression has been achieved if 
P < S' = AMAL; the "compression factor" is S'/P > 1. 

3) Form the PxS' "measurement matrix" M comprising those rows of the S' x S' matrix U^ (see (OTT l) 
whose indices correspond to the TF lag positions (m, /) chosen in Step 2. 

4) Compute an estimate r of r from a^^) by means of BP (l34l) . i.e., r = argmin]y[j./=a(^) Ik' 111 • 

5) From f, calculate Rx,cs['>T',k] = C{r}[n,k] according to (|36] |. This step can be implemented efficiently 
by two successive 2D FFT operations. 

Based on the error bound (l35T l. the compressive RS estimator Rxfisin, k] can be expected to be close to the 
noncompressive basic RS estimator RxM^^vi^jk] in (IT6l ) if the subsampled RS estimate i?x,Mvu[pAn, f/AA;] 
is approximately if -sparse and the index set Q in (1351 ) is chosen as described below (1351) . In the next section, 
we will derive a bound on the approximation error (MSE) that is formulated in terms of certain parameters 
depending on second-order statistics of the process X[n], including the RS, Rx[n,k]. 

We note that from an algorithmic viewpoint, our compressive RS estimator Rx,cs[n,k] is similar to the 
compressive TF representation proposed in [|54l . [|55l . However, the setting of Il54l . Il55l is that of deterministic 
TF signal analysis (improving the TF localization of the Wigner distribution), rather than spectral estimation 
for nonstationary random processes. 
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V. MSE BOUNDS 
In this section, we derive an upper bound on the MSE of the proposed compressive RS estimator Rx,cs [n,k], 

ecs = E{\\Rx,cs-Rx\\l} = Yl mRx,cs[n^k]-Rx[n,kf}, 

n,fce[iV] 

under the assumption that X[n] is a circularly symmetric complex Gaussian nonstationary process. We do not 
assume that the EAF Ax[m, I] is exactly supported on some periodic lag rectangle A = {— M, . . . , Mjjy x 
{-L, . . . , L}^ with < M < lN/2\ and < L < [iV/2j . 

A. Parameters 

Our MSE bound depends on three parameters of the second-order statistics of the process X[n], which 
will be defined first. 

• As a measure (in the broad sense) of the sparsity of Rx[n, k], we define the TF sparsity moment 

2 

(37) 



(w) A ^ 

a — 



"" WRxWl 



y^ w[n,k] \Rx[n,k]\ 

n,ke[N] 



where w[n,k] > is some suitably chosen weighting function and ||-Rx|l2 = SnfceWl l-^-'s^i"^' ^]| (i-^-' 
the norm is taken over one period of Rx[n, k]). In particular, for w[n, k] = 1, a)^ = ||^x||^/||^j!s:|l2- 

For another way to measure the TF sparsity, let us first denote by 

Rx,Myv[n,k] = E{Rx,MYv[n,k]} (38) 

the expectation of the basic RS estimator Rx,m^i][i^, k] in (fT6l ). It follows from ([8]l that RxM^vin, k] is 
a smoothed version of the RS, i.e., 

^x,Mvu[ra, A;] = — ^ ^MWv[n-n\k-k']E{Rx[n,k']} = — ^ ^M\Tj[n-n,k-k']Rx[n,k'] , 

n',k'e[N] n',k'&[N] 

(39) 
where E{Rx[n, k]^ = Rx[n, k] has been used in the last step. Due to (fTTI ). the smoothing kernel is given 

by 

^Mvu[n,fc] = ^ J] 0Mvu[m,/]e-^'t('=™-«O 

- ^ E lAm,l]e-^'-^^'--^^ (40) 

m,l&[N] 
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M L 

rn=—M m=—L 

Because of the smoothing, the number of significantly nonzero values of -Rx,mvu[^) ^] is generally larger 
than the number of significantly nonzero values of the RS Rxin, k]. However, for an underspread process, 
the RS is inherently smooth, which implies that the smoothed RS is close to the RS itself. Therefore, for 
an underspread process with a small number of significantly nonzero RS values, we can expect that also 
the smoothed RS consists of only a small number of significantly nonzero values. Let us denote by Q{K) 
the set of indices (p, q) G [AL]x[AM] of the K largest (in magnitude) values of the subsampled expected 



RS estimator, Rx,Mwi][p^n,qAk]. Let g{K) = [AL] x [AM] \ g{K), and note that \g{K)\ = S'-K. 



We then define the TF sparsity profiler 



|2i IP ^T2^n£>__ r„A^ „ALi|2i 



I^X||2 



with V - E{|(R)p+i,,+ir} = N^E{\Rx,MVv[P^n,qAk]\'} . (41) 
For later use, we note that 



J^V = E{||r^(^)||^} = E{||r-r^(^)||^}, (42) 

{p,q)£g{K) 



where r^(^) (r^(^)) denotes the vector that is obtained from r = vec{R} by zeroing all entries except 



the S' — K (K) entries whose indices correspond to the indiceslj {p,q) G G{K) i{p,q) G g{K)). 
The "TF correlation width" of X[n\ can be measured by the EAF moment IfTOl . ||24| 



,W A 



\\Axf 



Y^ ^[m,l]\Ax[m,l]\\ (43) 



2 m,l£[N] 



where tl;[m, I] is some weighting function that is generally zero or small at the origin (0, 0) and increases 

with increasing \m\ and |/|, and ||Aj>c||2 = ZlmiefM l^xi^T-i^ll = Il-^x|l2- For ^^ underspread process 

X[n] and a reasonable choice of the weighting function ^/J['m, I], rrr^ is small (^ 1). 

^We note that this definition is different from that in II56I . 

^For convenience, though with an abuse of notation, we denote by Q{K) both a set of indices k of (r)^ and the coiTesponding set 

of 2D indices (j>,q) of (R) , -^ , j = (unvec{r})p+i,5+i or equivalently of /Jx,Mvu[pAn, gAfc]. Thus, depending on the context, 

we will write k G G[K) or (p, q) G Q{K). 
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B. Bound on the MSE of the Basic RS Estimator 

Our bound on the MSE ecs = E|||^x,cs — ^^IL} ^^ ^ combination of a bound on the MSE of the basic 
(noncompressive) RS estimator Rx,M"^v[n, k] and a bound on the excess MSE introduced by the compression. 
First, we derive the bound on the MSE of the basic RS estimator, 

e = E|||i?x,Mvu — -Rxllaj • 
As in Section Unl we use the decomposition 

e = B^ + V, (44) 

with the squared bias term B^ = ||E{/?x,mvu} — ^x||2 ^"d the variance V = E{ ||i?x,Mvu — E{^x,Mvu}||2}■ 
ij Bias: An expression of the bias term is obtained by setting iplrrijl] = (j)M\i][m,l] = Ij([m,l] in (fT2l ): 

^' = E \ilA[m,l]-l)Ax[m,l]f = ^ I^[m,l]\Ax[m,l]\\ 

m,le[N] m,le[N] 

where I-^^im, I] is the indicator function of the complement A of the effective EAF support region A = 

{-M, ..., M}^ X {-L, ..., L}^, i.e., 



0, otherwise. 
We can write B^ in terms of the EAF moment (1431) with weighting function tplTn, I] = /^[jn-, /]: 

B^ = WAxWlm'i^^ = \\Rx\\lm'i^\ (45) 

Note that m-^-^ = 0, and thus B^ = 0, if the EAF Ax [tti, I] is exactly supported on A. 

2) Variance: In what follows, we will use the (scaled) discrete TF shift matrices 3jn,i of size NxN whose 
action on x € C^ is given by 

with {n)j^ = n mod A^. The family of TF shift matrices {3m,i},^ jgrjvi i^ considered in some detail in Appendix 
A. Using 3m,i, ^x,Mvu["^i ^] can be written as a quadratic form in x. In fact, starting from (IT6] l and using 
(1791) . we can develop Rx,M^i][iT'Tk] as follows: 



M L 

Rx,Myv['n,k] ^ ^ ^ y^ Ax[m,l]e 



O 1 V^ V^ ^ r_ 71 .-j^{km-nl) 



N 

m=-M l=-L 



17 



M L 

^ m=~Ml=-L 

/.ML \ 

= (-" ^ E E^^"^'"^-"^J-' 



Setting 



this becomes 



M L 

Cn,^-^ E E^'"^'™~"'^J-.^' (46) 



^x,Mvu[n,fc] = (xx^,C„,fe> = trjxx^C^J = x^C^^x. (47) 



Note that the matrix C„ ^ is not Hemiitian in general. 

Splitting iJx.Mvul'^, ^] into its real and imaginary parts, we have 

vaT:{Rx,MYu[n,k]} = var{^{Rx,Myv[n,k]}} + vaT:{Q{Rx,MYv[n,k]}} . (48) 

The real part of -Rx,mvu["^) k] can be developed as 

3f?{i?x,Mvu[n, A;]} = ^ (i2x,Mvu[n, A;] + i?3f,Mvu[^> k]) 1* ^ (x^C^feX + x^C„,fcx) = x^ci'^^x , (49) 

with the Hermitian matrix 

CS = I «k + Cn,k) . (50) 

Similarly, we obtain for the imaginary part 

9{i?x,Mvu[n,A;]} = x^C^'^^x, (51) 



with the Hermitian matrix 



^n,k - ^ (Cn,fc - Cn,k) ■ (52) 



Inserting ( |49l ) and (ISTI ) into (1481 ) and using a standard result for the variance of a Hermitian form of a circularly 
symmetric complex Gaussian random vector ll57l . we obtain 

var{i?x,Mvu[n, k]} = trjC^rx^irx} + tr{C« rxC« Fx} • (53) 

Using this expression, we next derive an upper bound on y = Ej ||i?x,Mvu — Eji^x.MvullL}- ^^ have 

y= ^ E{|^x,Mvu['n-, /c]-E{^x,Mvu["-, A:]}| } 

n,k£[N] 
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= ^ vaT{Rx,MYv[n,k]} 

n,k&[N] 

^ E trlcSr.cSr,} + Y. tr{c« rxc« r-} • (54) 

n,ke[N] n,ke[N] 

It is then shown in Appendix B that 

V = Yl \Ax[^,l]\\[^,l], (55) 

m,l&[N] 

with 

1 ML 

X[m,l] = j^ E lA[m',l']e^'^'~^^'-^'^ = j^ ^ ^ ^^ " ^'"'""'^ ' (^6) 

m',/'e[Af] m'=-Ml'=~L 

We can bound the magnitude of x[m, /] according to 

ML „ 

|XK/]| <^ E E le^'-^''"'-™'')] = -(2M + 1)(2L + 1) = -. 

m'=-A/P=-L 

Combining with (1551 ) leads to the following bound on V: 

V< Y \Ax[^J]\'\x[m,l]\<^ Y |^x[m,/]|'i|px|l2. (57) 

m,l&[N] m.,l&[N] 

3) MSE: Finally, the desired bound on the MSB e = E{ H-Rx.mvu — ^^Ho} ^^ obtained by inserting (1451 ) 
and (|57] ) into the expansion (l44b : 



,2^„ . «i...,,^.j«.|p,ii= = p,ii:^/,„to) + | 



s = B^ + V < \\Rx\\im\i--> + -\\Rxl\i = II-Ra-112 mF'+T7r 08) 



This bound is small if X[n] is underspread, i.e., if m-^-^ <^ 1 and S -^N. 



C. Bound on the Excess MSE Due to Compression 

The excess MSE caused by the compression is given by 



Ae = E{||i?x,cs — -Rx,Mvu||2} 



Because of the Fourier transform relations (l32i and (136) . we have 



Ae = lE{||r-r||^}. (59) 
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As in Section IIV-C[ let K denote a nominal sparsity degree that is chosen according to our intuition about the 
approximate sparsity of i?x,Mvu[pAn, gA/c] and, equivalently, r. We assume that the number P of randomly 
selected AF samples is sufficiently large so that (1351 ) is satisfied, i.e.. 



12 - ^11 111' 



(60) 

for any index set Q of size \Q\ = K. (A sufficient condition is (1331).) An intuitively reasonable choice of K 
and Q can be based on the smoothed RS Rx,MNv[n, k] = E|/?x.mvu["^5 ^]} in (l38l ). (l39l ): we choose K as the 
number of significantly nonzero values Rx,MW\]\p^'n, gAfc], and Q = Q{K) of size K as the set of those indices 
of r such that the corresponding values i?x,MvubAn,gAA;] are the K largest (in magnitude) values. Thus, 
yG{K) comprises those K values Rx,MNv[p^'n'-,q^k] for which the corresponding values Rx,MNv[p^'n,q/^k\ 
are largest (in magnitude). 

Based on this choice, we will now derive an approximate upper bound on the excess MSB Ae. Inserting 
(l60l ) into ^9\, we obtain 



D 



2 



A6< -^E{||r-r5(^)||^}. (61) 

Using the inequalitjzl || • H^ < || • \\q\\ ■ \\l, we have ||r-r^(^)||5 < ||r-r^(^)||J|r-r^W||2 < {S'-K)\\y- 
r^(^)|| , and thus (|6TI ) becomes further 
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-l2 ^ rr^ fQl l^\ n2 



(5^-K)^ II a(;,)||2, m {S'-K)D^ ^ m {S'-K)D^ . 

¥k ^lF-^ II2I - gf^ 2^ V - s^ ll^xlls f^x(A). (62) 

{p,g)egM 

In what follows, we will derive an approximate expression of hp^q = E{|(R)^^ ,^| } in terms of 
Rxin, k\; this expression will show under which condition ax{K) oc ^, \g q(k) ^p-i ^^ small. We have 

hp,q = E{|(R)p+i,g+il } 

= vaT{mnp-,i,,+i}}+MHi^)p+i,,+i}} + |E{(RUi,,+i}|'. (63) 

Using (l27l) and ( 1791 ). we can express (R)^^ _^_i as a quadratic form: 



M 

(R 



^ E E^^t"^'^]^"''"^^"^^ 



p+l,g+l — 

m=-M I 



^Indeed, the £1 norm of an arbitrary vector z can be expressed as |[z||-^ = z^a(z), where a(z) is given elementwise by (a(z))(. = 
Zk/\zk\ for 2fc / and (a(z)), ^ for Zk = 0. Clearly, |[a(z)||^ = |jz||o, and thus ||z|[J = (z^a(z))2 < \\z\\l ||a(z)||^ = ||z|[^ \\z\\g, 
where the Cauchy-Schwarz inequality has been used. 
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M L 

= ^E E (--''' J-.') ^"'"^^ 

m=-M l=-L 

= (xx^,Tp,,) 
= trjxx^T^J 



x^T^x 



9^ _ P^ 
AM AL y 



with 



Af L 



%. = ^E E^^"^^ 



AM Ai / T , 



m=-M l=-L 

Note that the matrix Tp g is not Hermitian in general. Inserting (|64] | into (l63T l then yields 

|2 



/l. 



'P.? 



varjx^TWx} +var{x^T«x} + |E{x^Tp^,x}|^ 



(64) 



(65) 



with the Hermitian matrices 



rp(R) A l/rpi? , rp \ rp(l) A J_ /rp// _ rp \ 



(66) 



Using standard results for the variance and mean of a Hermitian form of a circularly symmetric complex 
Gaussian vector ll57l . we obtain further 



V = trjTWrxTWrx} + tr{TWrxT«rx} + |tr{rx<j|^ 



(67) 



There does not seem to exist a simple closed-form expression of (l67l ) in terms of the EAF Ax[m,l] or 
the RS Rxin, k]. However, under the assumption that the process X[n\ is underspread and the effective EAF 
dimensions L, M (cf. (I36l )) are accordingly chosen to be small, the following approximation is derived in 
Appendix C: 



El _ I 2 



y Rx[n,k\'^MVv[n—p^'n,k — q/S.k] 

n,k&[N] 



, (68) 



n,fce[Af] 

where, as before. An = N/{AL) and Ak = N/{AM). Comparing with (|39l ) and noting that $mvu[— "^, —k] = 
^Mvui'^, k], it is seen that the second term on the right hand side of (l68T l is N"^ \Rx,Mvv[p^n, qAk]\ . Using 
the inequality ||-||2 < Irlli llSllI to bound the first term on the right-hand side of (I68l l. and using a trivial 
upper bound on the second term, we obtain 



hp,i ^ N 



y \Rx[n,k]^M\u[n—pAn,k — qAk]\ 

n,k&[N] 



+ 



E \Rx[n,k]^MYv[n-pAn,k-qAk]\ 

n,k&[N] 
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(iV+1) 



n,ke[N] 



(69) 



Here, X]n ke\N]\^x[n, k] $mvu[^— pAn, A; — gAA;]| can be interpreted as a local average of the RS modulus 
|^x[^;^]| about the TF point {pAn,qAk). Thus, the (approximate) upper bound (|69] | shows that /ipg is 
small if Rx[n, k] is small within a neighborhood of (pAn, (/A/c) or, said differently, if {pAn, qAk) is located 
outside a broadened version of the effective support of Rx[n, k]. The broadening is stronger for a larger spread 
of ^Mwv[n,k]. According to (|40] |. <^mvu[^5^] is the 2D DFT of the indicator function Ij^mJ], and thus the 
broadening depends on the size of the effective EAF support A; it will be stronger if A is smaller, i.e., if the 
process X[n] is more underspread. Since a stronger broadening implies a poorer sparsity, this demonstrates 
an intrinsic tradeoff between the underspreadness and TF sparsity of X[n]: better underspreadness implies a 
smaller effective EAF support A, whereas better TF sparsity requires a larger A. 



With this "broadening" interpretation in mind, we reconsider dx{K) oc Yl,( 



■hria in the bound 



(l62l ). Recall that G{K) was defined as the set of those indices of r such that the corresponding values 
i?x,Mvu[pAn, qAk] are the K largest (in magnitude). Therefore, a small ax{K) requires that K is chosen such 
that KAnAk is approximately equal to the area of the broadened effective support of Rxin, k], because then 

0. 



J2nke\N]\^x[n,k] '^MYv[n-pAn,k-qAk]\ ^ for {p,q) £ g{K) and thus dx{K) oc Y. 



{S'-K)D^ 



Using ( 1691 ). we can upper-bound the MSB bound in (1621 ). Ae < '■'^ s^K^ Sf )eG(K) '"P'I 
in a simpler (but generally looser) upper bound. Indeed, we have 



(p,q)egiK)^P,<i 

hp q, which results 



{p,q)&g{K) 



< 



ip,q)egiK) 



y \Rx[n,k]^uwv['n—pAn,k — qAk]\ 

nMW] 



< (N+1) 



(AT + l) 

(iV+1) 



{p,q)&g{K) 



y \Rx[n,k]^MWv['>T'—pAn,k — qAk]\ 

n,k&[N] 



{N+1) 



y^ y^ \Rx[n,k]\\^M\i][n-pAn,k-qAk]\ 

(p,q)e'g(K)'nMlN] 

Y^ \Rx[n,k]\ Y^ \<pMMv[n-p^n,k-qAk]\ 

n,ke[N] {p,q)eg(K) 

-| 2 

y \Rx[Ti,k]\w^[n,k] 

n,k&[N] 



(70) 



where 



< II • 11^ was used in the step labeled with (*) and 

w^[n,k] = y^ \^M\\]['>T'—pAn,k — qAk]\ 
ip,g)^G{K) 



(71) 
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Comparing with the definition of the TF sparsity moment aj^ in (l37l ). it is seen that the approximate bound 

(TTOl) can be written as 

J^V, < {N+l)\\Rx\\la^r^. (72) 

(p,g)egM 

Inserting (1721 ) into (l62l ) then gives the approximate MSB bound 

Ae ^ ^^^^^§^(^^+1) Pxll^a?^^ (73) 

A small excess MSB Ae can be achieved if the TF sparsity moment cj^* oc ^„ ke\N] \^x[n, k]\ w^[n, k] is 
small. This, in turn, is the case if the RS Rx [n, k] is negligible within the effective support of the TF weighting 
function w^[n, k]. Due to (TTT]) . the size of the effective support of w^[n, k], which is concentrated around the 
points |(pAn, qAk)^ , )g c(k) ' ^^ ^^^ larger than S'—K times the size of the effective support of <1>mvu[^, k] 
(recall that \GiK)\ = S' -K). Because of the DFT expression (gOll and the fact that \[Nf' r\ A\ = S (see 
(l22l) ). the size of the effective support of $mvu['^> k] within one period [N] can be estimated by N'^/S. Thus, 
for a small TF sparsity moment o^* , the RS Rxin, k] should effectively vanish on a region of size at least 
{S'-K)N^/S > {S- K)N^/S = N^ - KN^/S. Since typically S' ^ S, implying that {S'-K)N^/S « 
N'^ — KN'^/S, it follows that the size of the effective support of the RS Rxi'n, k] should not be larger than 
KN'^/S. Note that K was defined as our prior intuition about the number of significantly nonzero values 
i?x,Mvu[pAn, gA/c]; furthermore, N'^/S is related to the TF undersampling in i?x,Mvu[pAn, gAA;] because 
(for S' K, S) it is approximately the ratio of the number of samples {i?x,Mvu["^; ^]} h(z\N] (which is N"^) to 
the number of samples {-Rx,Mvu[pAn, gAA;]}pgj^^j^gj^^_^j (which is S"). 

D. Combining the Two MSE Bounds 

We will now combine the bound (1581 ) on e = E{||/?x,mvu — ^^IL} ^^^ ^^^ bound (l62l ) or (1731 ) on 
Ae = E|||i?x,cs — -Rx.MvulL} ii^to a bound on the MSB ecs = E|||i?x,cs — ^^Ho} °^ ^^^ proposed 
compressive RS estimator Rx,cs[n, k]. To this end, let us define the norm of a random process Y[n, k] that is 
A^-periodic in n and k as 



riiR = V^iii^ii'} = J E mY[n,kf}. 

\ n,k&[N] 

For two such random processes Yi[n, k] and Y2[n, k], the triangle inequality lISTI states that 

i|yi + ^2iiR<riiiR + r2iiR. m 

The estimation error of the compressive RS estimator can be expanded as 
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Rx,cs [n, k] - Rx [n, k] = Rx,cs [n,k]- -Rx,mvu ["-, k] + -Rx,mvu [n, k] - Rx[n, k] 

= Yi[n,k]+Y2[n,k], 

where we have set Yi[n,A;] = Rx, Mwin, k] — Rx[n, k] and 1^2 [^,^1 — Rx,cs[i^, k] — Rx,Mwv[n, k]. Hence, the 
MSE of the compressive RS estimator can be rewritten as 

ecs = E{\\Rx,cs-Rx\\l} = E{||yi + y2|l2} = \\Yi + Y2\\l. 
Using the triangle inequality (1741 ). we then obtain the bound 

ecs < (ri||R + r2||R)'. (75) 



Recognizing that ||1i||r = y E{||i?x,Mvu - -RXII2} = \/e and ||12||r = yEJUi^x 
the bound (1751) becomes ecs < (\/e + "\/Ae) and, finally, 

ecs < (^ + \/Ai)^ 



Rx. 



MVU 



^, 



Inserting the bounds (l58l ) on e and (l62l) on Ae then results in the following overall bound on ecs^ 

2 



ecs < \\Rx\ 



(/_) S {S'-K)D^ ^ ,^^. 



Alternatively, using the approximate bound (1731 ) on As instead of (l62l ). we obtain the simpler (but looser) 
approximate bound 

i2 



ecs 



< 



\Rx\ 



'- +iV + 



^^^^l^J^-^^^iN^l)^?^^ 



We note that our bounds on Ae are based on the CS bound (lT9l ) together with (ITSl ). which is known to be 
very loose. Thus, for a given nominal sparsity degree K and a given number of measurements P satisfying 
(l33l) (cf. (fTSl)). our upper bounds on Ae and, in turn, on ecs will generally be quite pessimistic, i.e., too high. 
However, the bounds are still valuable theoretically in the sense of an asymptotic analysis, because they show 
how the MSE decreases with increasing underspreadness (expressed by a smaller moment m)^^ ) and with 
increasing TF sparsity (expressed by a smaller moment cr^* ). 

VI. Numerical Study 

We will assess the performance of our compressive spectral estimator in a setting that is inspired by a 
cognitive radio application. In a cognitive radio system, a given transmitter/receiver node has to monitor a large 
overall frequency band and determine the unoccupied bands that it can use for its own transmission ||3]-||5]. 
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A. Simulation Setup 

We consider a single active transmitter employing orthogonal frequency division multiplexing (OFDM) ||58]| . 
||59l . which is a modulation scheme used, e.g., for wireless local area networks (WLAN) ||59l . Il60l . digital 
video broadcasting (DVB) 11611 - 11631 . and long term evolution (LTE) cellular systems | l64l . In our simulation, 
the ODFM transmission uses L = 64 subcarriers and a cyclic prefix whose length is 1/8 of the symbol length. 
Each subcarrier i G [L] transmits a symbol Sj that is randomly selected from a QPSK constellation with 
normalized symbol energy £'s = |sip = 1. All QPSK symbols are equally hkely, and the different subcarrier 
symbols Si are statistically independent. The OFDM modulator uses an inverse DFT of length L = 64: to map 
the frequency-domain transmit symbols Si into the (discrete) time domain; this is followed by insertion of 
a cyclic prefix. Assuming an idealized, noise-free channel for simplicity, the resulting transmit signal is also 
observed by the receiver. However, we assume that our receiver monitors an overall bandwidth that is twice the 
nominal OFDM bandwidth, B. This corresponds to a twofold oversampling, i.e., a sampling period of 1/{2B), 
and can be easily realized by using an inverse DFT of length A'^, = 2L = 128. The lengths of an OFDM 
symbol and of the cyclic prefix are then given by Nf^ = 128 and A^^^p = 128/8 = 16 samples, respectively. To 
keep the simulation complexity low, we assume that a single OFDM symbol is transmitted, with silent periods 
before and afterwards. Thus, the received time-domain signal (discrete-time baseband representation) is given 
by 



X[n] 



ie[L] 
0, 



else. 



Here, no denotes an arbitrary but fixed time offset. In our simulation, we used uq = Ncp and considered X[n] 
for ne [N] with A^ = 512. 

Because of the random Sj, X[n] is a nonstationary random process. The RS and EAF of X[n] are easily 
obtained from, respectively, (S) and ([1} as 



Rx[n,k] = < 



ie[L] 
0, 



V ■ N 2L 



else; 



Ax[m,l] = < 



ie[L] 

-j^l^ A*. 



I 



m, 



-cp '-, ^ 

0, 



— p-^^ir™ 



e-^-™, me{-N,p-N, + l,...,0} 



'N 



else, 
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Fig. 1. TF representation of the OFDM process X[n]: (a) RS Rx[n, k], displayed for (n, k) e [N] x {-iV/2, . . . , N/2 - 1}, with 
iV=512; (b) EAF Ax[m,l], displayed for {m,l) G {-iV/2, . . . , iV/2 - 1}^ 



^n-1 



where dir(n, ^) = ^^'=0^''^^'^' ~ e^^^^^ ^^ '^sin^re\ ■ Note that the expression for ^xi^T-,^] requires that 
Ncp + N, < N/2, a condition that is fulfilled in our simulation since 128 + 16 < 512/2. The RS and EAF are 
shown in Fig. [T] From this figure, we can conclude that the process X[n\ is reasonably TF sparse but only 
moderately underspread (the latter observation follows from the fact that Rx [n, k] is not very smooth). Note 
that the TF sparsity could be further improved if we considered longer silent periods before and/or after the 
OFDM symbol, and if we considered a wider band (i.e., if we used an oversampling factor larger than 2). 
For the design of the compressive RS estimator Rxfisi^, k] in (l36l ). we used M = 3, L = 7 and AM = 8, 



AL = 16. This corresponds to choosing the effective EAF support (see Q) as ^ = {—3,..., 3} 



512 



{—7, . . . , 7}5]^2' of si^^ '^ — (2^^ + 1)(2^ + 1) = 105; furthermore, the size of the extended effective EAF 
support A' is S' = AMAL = 128. For an assessment of the TF sparsity of X[n], we consider /ipg = 
A^^ E| |i?x,Mvu[pAn, gAA;]| }, which underlies the TF sparsity profile axiK) in (|4T] ). Let {p,q)j. with r G 
{1, . . . , S'} be the TF index of the rth largest (in magnitude) value of the set {i?x,Mvu[pAn, gA/c]} , iprAn [amI' 
where, as before, Rx,M\\j[n,k] = E{Rx,MV\j[n,k]} = j^Y.n',k'&[N]'^M\\j[n-n',k-k'] Rx[n',k'] (see 
(l39l)). In Fig. [2j we show the values /i(„ „■) along with the corresponding approximations 



-here denoted 



h 



{p,q) — ^s a function of the index r. It is seen that h(^pq-j is close to zero for r larger than 15. Furthermore, we 
can conclude that the ordering of the values i?x,MvubAn, gAA;] according to decreasing magnitude matches 
the ordering of the values hp^g very well. Thus, for TF positions {pAn,qAk) for which |i?x,MvubAn, gAA;]| 
is large, we can expect that also hpg is large. Finally, it is seen that the curves representing hrp^g\ and /i(p^q) 
coincide, which shows that the approximation (l68T l is very accurate. 
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Fig. 2. h(^p_q-) normalized by /i(p,5) and the corresponding normalized approximation /i{p,g) according to ( I68t versus r. 
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Fig. 3. Averages and single realizations of RS estimators: (a) RS of the process X[n], (b)-(d) averages and (e)-(g) single realizations 
of the RS estimates obtained with the following estimators: (b), (e) noncompressive estimator Rx,u'^v[n,k] (compression factor 
S' / P = 1); (c), (f) compressive estimator Rx,cs,[n;k] with S' / P — 2; and (d), (g) compressive estimator with S'/P ^ 5. All TF 
functions are shown for (n, k) £ {-150, . . . , 361} x {-N/2, ..., N/2 - 1}, with A^ = 512. 



B. Simulation Results 

We now consider the estimation of the RS Rx[n,k] from a single realization of X[n] that is observed 
for n G [512]. To evaluate the estimation performance, we generated 1000 realizations of the QPSK symbols 
Si, i^ [64] and computed the corresponding realizations of X[n]. In Fig. [3l we show single realizations along 
with the average of the compressive RS estimates Rx,cs[''T', k] obtained for the 1000 realizations of X[n\, for 
compression factors S'/P = 1, 2, and approximately 5 or, equivalently, P = 128, 64, and 25 randomly located 
AF measurements. The true RS is also re-displayed for easy comparison with the average estimates. The case 
S'/P = l corresponds to the basic RS estimator i?x,Mvu['i^5 ^] in (IT6l ). We see that already in this case, there 
are noticeable deviations from the true RS. In fact, the average of the 1000 basic RS estimates Rx,M^i][n, k] 
closely approximates the expected basic RS estimator i?x,Mvu[";^] = E{i?x,Mvu["';^]}. which according to 
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Fig. 4. Empirical normalized MSB, squared bias, and variance of the compressive RS estimator -Rx,cs[?^, fc] versus the compression 
factor S' IP. 



( 1391) is a smoothed version of the RS. This smoothing leads to a significant deviation from the RS, because 
the RS itself is not very smooth. The limited smoothness of the RS corresponds to the fact that the process 
X\p\ is only moderately underspread. For compression factor S' jP = 2, there is no visible degradation of the 
average estimate relative to the basic estimator For S' jP ss 5, a small degradation is visible. 

For a quantitative analysis of the degradation caused by the compression, we show in Fig. |4]the empirical 
normaUzed MSE (NMSE) of the compressive estimator Rx,c^\n^ ^] versus the compression factor S' / P. The 
NMSE is an empirical, normalized version of the MSE ecs = Ej ||^x,cs— ^xjL}' with the expectation replaced 

II ~ I|2 

by the sample average over the 1000 process realizations and with normalization by ||i?x||2- Iii the same figure, 
we also show the empirical normalized versions of the squared bias term B^^ = ||E{i2x,cs} "^^IL ^"^^ ^^^ 
variance Vcs = E|||i?x,cs— EliJx.csllL}' again with normalization by H^xIL- (Recall that ecs = ^cs+^s-) 
These results demonstrate a "graceful degradation" with increasing compression factor S' / P. Again, the result 
for S' /P = 1 corresponds to the basic RS estimator i?x,Mvu[^5 k]- We did not plot the MSE bounds derived 
in Section FVl because they are much larger than the empirical MSE. As mentioned in Section IV-Di this lack 
of tightness is mostly due to the notoriously loose Il53l CS error bound used in (|60] | (combined with (|33]|). 

VII. Conclusion 

For estimating a time-dependent spectrum of a nonstationary random process, long-term averaging cannot 
be used as this would smear out the time-dependence of the spectrum. However, if the spectrum as a function of 
time and frequency is sufficiently smooth, which amounts to an underspread assumption, a local TF smoothing 
can be used. In particular, the RS of an underspread nonstationary process can be estimated by a local smoothing 
of a TF distribution known as the Rihaczek distribution. 
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In this paper, we have considered the practically relevant case of underspread processes that are approx- 
imately TF sparse in the sense that only a moderate percentage of the RS values are significantly nonzero. 
For such processes, we have proposed a "compressive" RS estimator that exploits the TF sparsity structure 
for a significant reduction of the measurements required for good estimation performance. The measurements 
are values of the discrete AF of the observed signal at randomly chosen time lag/frequency lag positions. 
Our overall approach is advantageous if dedicated hardware units for computing values of the discrete AF 
are employed, and/or if the AF values have to be transmitted over low-rate links or stored in a memory. 
The proposed compressive RS estimator extends a conventional RS estimator for underspread processes (a 
smoothed Rihaczek distribution using a minimum variance unbiased design of the smoothing function) by a 
CS compression-reconstruction technique. For the latter, we used the BP because it is supported by a convenient 
performance guarantee (a bound on the £2 norm of the reconstruction error); however, other CS techniques 
can be used as well. 

We provided upper bounds on the MSE of both the minimum variance unbiased RS estimator and its 
compressive extension. The MSE bound for the compressive estimator is based on the error bound of the BP, 
which is known to be quite loose. Therefore, the MSE bound for the compressive estimator is usually quite 
pessimistic. However, the MSE bound is still useful theoretically, since it reveals the asymptotic dependence 
of the estimation accuracy on the underspreadness and TF sparsity properties of the process. Numerical 
experiments showed that for a typical scenario, our compressive estimator works well up to compression 
factors of about 5. 

We considered the RS because in the discrete setting used, the RS is the simplest time-dependent spectrum 
from a computational viewpoint. However, for underspread processes, the RS is very close to other important 
time-dependent spectra such as the Wigner-Ville spectrum and the evolutionary spectrum. Therefore, the 
proposed RS estimator can also be used for estimating other time-dependent spectra if the process is sufficiently 
underspread. Finally, the proposed RS estimator can also be used for estimating the EAF, based on the 2D 
DFT relation connecting the RS and the EAF. 

Appendix A: TF Shift Matrices 

We consider the family of (scaled) discrete TF shift matrices {"^ m,i} ^^ ipij^] of size NxN whose action 
on X G C^ is given by 

(Jm,/x)„+i = -^ (x)(„_„)^ + i e^'i^'", ne[N], (76) 

with {n)j^ = n modiV. These matrices can be written im,i = -7^ M^Tm, where M; is the diagonal NxN 
matrix with diagonal elements 1, e-''""', . . . ^e^^''^^^^' and T^ is the circulant NxN matrix whose entries 
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(Tm)„ n' ^^ given by 1 if {n — n')^ = (m)^ and otherwise. It can be easily verified that the set {Jm,^}^ iptj^] 
forms an orthonormal basis for the linear space of matrices C^^^ equipped with inner product (A,B) = 
tr{AB^}, i.e., 

(Jm,«, Jm',p) = S[m-m]j^6[l-l']j^ (77) 

and 



A= Y, {A,Jm,l)J 
m,l(^[N] 



ml 



for all A G C 



NxN 



It can furthermore be shown that the EAF in ([1} and the AF in (l3j can be written as 



Ax[m,l] = VAf(^x,J„^,^> 
Ax[m,l] = ViV(xx-f^,J„,,), 

where Tx = E{xx^}. Thus, according to (iTST l. we have the expansions 



XX 



H 



m,«e[Af] 
m,iG[Af] 



Finally, from (1761 1. one can deduce the following relations: 

1 , 



"ml "m'l' 



N 



i+m'A + Ve ' « 






and, in turn. 



''ml — -'-m-l K : 



''n,k"m.,l"n,k — t^t "m,,l^ 



(78) 



(79) 



(80) 



(81) 
(82) 

(83) 



Appendix B: Derivation of Expressions ([551) and (l56l ) 

We will derive (ISST i and (l56l ) from (l54l l. Our derivation will be based on expansions of C|^ l and C^^ into 
the TF shift matrices 3m,i- Using (l50l ). (I46l ). and (l82l ). we have 

ML ML 



n,k 2ViV 



'mJ 



(82) 1 



2VN 



--Ml=-L m=-Ml=-L 

ML ML 



Lm=-Af i=-L 



m=-Ml=-L 



with 



1 



2ViV 



M L 



M L 



Y^ Y^e^'-^ikm-nDj^^^^-j^-^ml ^ Y^ ^ ^^^-^ikm-nl) j 



ml 



m=-Ml=-L 
M L 



m=-M l=-L 



^ E E e^' " ^'™-"'^ (e-^'^'"' + 1) Jm,i 



2VN 



m=-M l=~L 



m,l&[N] 



.(R) A 
"m,l 



1 



2VN 



X^[m,/](e-^'^'"' + l) 



In a similar way, we obtain from (|52] | 
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(84) 



(85) 



.(I) 



, — > e-' « 

n,fe / y 



J^{km-nl) m 



c'1'3. 



I "ml 1 



with 



2jy/N 
The first term in (|54] | can then be written as 



c«, 4 ^:^x^[^,;](e-^-t-'-i). 



(86) 



(87) 



with 



E t^cSrx^fc} = tr{Dr^} = tr{Dr#} 



n,k&[N] 



D = E cSr.<) 

n,fce[Af] 
n,fee[Af] 



<11 Y^ V^ Y^ ^(R)^(R)*^j^[fe(m-m')-n(Z-r 

n,fce[Af] m,l&[N] m',l'&[N] 



iH 



/ . / , ^m,l'^m\V^ -im,li-X'im',l' 



E E 

^,l&[N]m',l'e[N] 



(R) (R)* 

m,l m',l' 



\ QJ N 

n,ke[N] 



j^[k(m-m')-n{l-l')] 



Jm,l^xJrn',l' 



N^5[m-m']j^5[l-l'] 



N^Y 



N I -I JV 



.(R)|2 



Jm. ;rvJ 



H 



^^l\ -im,li-X-im,l 



m,l&[N] 



in) Ar /^ Y^ Y^ I (R) |2 ;f r / ;/! 7 t t^ 



m,lelN]m',l'elN] 



(83) 



^/^ E E IcSRxKnJmVe 



-j^i-ml'-lm') 



m,l&[N]m',l'&lN] 



(88) 



(89) 
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and 



.ii- (p 1 



^x 



m,l&[N] 



(90) 



Inserting ^^ and (|90ll into ([88]) then yields 

E trKlr.c^rx} 



tr 



\ m,Lm',l',m",l"G\N] ) 



^",1"&[N] 

.(R)|2 



m,Z,m',J',m,",i"£[A''] 



(J„,,,,J„.,,,">'S'5[m'-m"]„<5[«'-P'] 



W.l^l^^fm'/'lPp-i^M'-''"') 






m,l,m',l'&[N] 



(91) 



(92) 



m,Z,m',Pe[Ar] 

where the relation 

Ax[-m,-l] = A*x[m,l]e-^^'^^ 

has been used in the last step. 

In a similar manner, using (|86] |. we obtain for the second term in (l54l) 



E tr{C« r.C« r.} = 5: |cS/ l^xK if e^-^-^--'^ . (93) 

n,fee[Af] m,l,m',re[N] 



Inserting ([93 and (|93]l into ^5^ then gives (|55] ): 



with 



m,Z,m',/'e[Af] m,Z,mye[Af] 

= Yl \Ax[m',ltx[m\l'], 

m',l'&[N] 



m'J'e[Af] 



Using ([851) and ([87]), we have 



1 



\ in,l\ ' \m,l\ i\j -^L 'J 



j;^^A[m,l] 



e~J¥"^i + 1 



+ 



-j^-^i _ 1 



2j 



vr 



iV 



cos I -T7^ + sin I —-:ml 



vr 



A^ 



1 



lj\^[m,l], 



(94) 
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and hence ( |94l ) becomes 



X[m,^ = j;^ "^ XjXm',l']e 



M 



J^im'l~l'in) 



m',l'&[N] 



N ^^ 

m'=~M l'=-L 



y^ f,i^ira'l-Vm) _ 



which is (l56ll. 



Appendix C: Derivation of Expression 



We will derive ^ from (l67T i. 

1) Expansions of TpJ and Tp/q: Our derivation will be based on the underspread assumption and on 
expansions of TpJ and TpJ into the TF shift matrices 3m,i- Inserting (1651 ) into the definition of TpJ in (l66l ) 
yields 



rp(R) 

p,q 



-i2-(i^-fi. 



~ 2 



M L 

E E 

rra=-A/Z=-L 
ML 

m=-M/=-L 



M L 



<i+ E E^^'"^^""^J"^>^ 



m=-M l=-L 



M L 



- 7 2^ mZ 



gTTi _ pi 



m=-M l=-L 



M L 



m=-Ml=-L 



a=-Ml=-L 1 

m=-A/Z=-L 

E,-27r('-22L_JiL) ,(R) _ 

m,le[N] 



with 



.(R) 



X^[m,/](e-^«"^' + l). 



In a similar manner, we obtain the expansion 



(95) 



(96) 



m,l(^[N] 



with 



i«, = ^X^[m,Z](e-^-^-'-l). 



ml 



2i 



(97) 



(98) 



For an underspread process X[n], the effective EAF support A = {— M, . . . , M}^ x {— L, . . . , L}jy is a 
small region about the origin of the (m, /) plane (plus its periodically continued replicas, which are irrelevant 
to our argument and will hence be disregarded). Looking at the expressions of t^i and tj^^ in ( |96l ) and ( |98l l. 
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we can then conclude from the presence of the factor Iyi[m,l] that t^^ and t^^ can be nonzero only for 
|?TT,/| ^ N. This means that in ( |96l ) and ( |98l ). we can approximate e"-' «"™' by 1, yielding 

C ~ ^lA[m,l] (99) 

4!/ « • (100) 
Using (ITOOl) in (|97]l yields 

T« « , (101) 



and thus (I67b approximately simplifies to 

V ~ trJTWrxTWrx} + |tr{rxTp^j|'. (102) 

2j First term in (I102I ).- We will now develop the two terms on the right-hand side of (1102b . The first term 
can be written as 

tr{Ti5rxTWrx} = (T^rx, (TWr^)^) . (i03) 

In order to find an approximation for this inner product, we use the following general result for the product 
C = AB of two NxN matrices A and B. The matrices A, B, and C can be expanded into the orthonormal 
basis Um,i}m,ie[N]' ■*!* respective expansion coefficients am,h hm,h and c^,/, e.g., A = Y.m,i&{N] 0'm,i^m,i- 
Then the c^,/ are related to the a^^i and 6^. / by the "twisted convolution" 1241 . Il33l . Il65l . Il66l 



Cm,i = -^^ V am',vhm-m;i-i'e-^'^'^'^^-^'y (104) 

This expression can be verified using (ISTT l. Let us apply it to the matrix product Tp^^ Fx- We have the expansion 

The expansion coefficients of Tp^^ are e-^^'^v^"^^ t\^\ (see (|95])); those of Ty are -4= Ax[m,l] (see (l80l)). 
Inserting these expressions into (I1041 l yields 









Ax [m — m', I — I' 



N 



g-j^m'{l-l') 



^ Y^ e^'2^(Si^-^)/4K/']^x[m-m',/-r]e-^«"*'('-''). (106) 



N 

m',l'e[N] 

where the approximate expression (|99] l was used in the last step. For an underspread process, because of the 
support of lA[fn, I] and the effective support of Ax[m, I], the terms in the sum (IIO6I 1 are significantly nonzero 
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only for \m'{l — l')\ ^ N. We can thus use the approximation e ^ «"'(' ') « 1 in (11061 ). which yields 

dp,,.^,;«-^ J2 lA[m',l']Ax[m-m',l-l']eMTKi-^) . (107) 

Next, we consider 

/rp(R)T-, \H (Ml V^ ,* ^H 

\-'-p,q^X) / J "'p,q-m,l''m,l 

m,l(^[N] 

~ Z-/ "'p,q;m,l''-m-l^ 

m,l&[N] 

= E 4,.;-m,-/J-,;e-^'^™'. (108) 

rri,Ze[Af] 

For an underspread process, again because of the support of I_A[m, I] and the effective support of Ax[m, I], it 
follows from (llOVI l that the coefficients dp^q-m,i are significantly nonzero only for \7nl\ <^ N. Hence, we can 
set e"-'^""' Ri 1 in (11081 ). which gives 

{T^^xf^ E dlg..m,-lJm,l. (109) 

m,l(^[N] 



In a similar way, we obtain from (I92l i the following approximation: 

^x[-m, -?] w A3^[m,Z], for|m/|<iV. (110) 

We now insert (llOSI l and (|1091 l into (11031) . and obtain 

^^\^p,q^X i-p^q^X J ~ ( / ^ dp^q-rn,l 'i m,l , / ^ "p,q;-m,',-i' Jm,',/' ) 

\m,«e[Af] m',l'£[N] I 

/ ^ "'p,q;m,l ^p,q; — m, — l ■ \^^^) 

m,lelN] 

From the underspread approximations (1107) and (IllOt . it readily follows that dp^q-^rn-i ~ d* -mi- Indeed, 
dp,q,.m,-i^ -^ E /^Kn^x[-m-m',-/-ne^-2'^(S^-^) 



fTTol 1 



ML 

^ ^^ m'=-Ml'=-L 
.ML / , ,\ 

^^ m'=-Ml'=-L 
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(*) 1 



7N ^ '^■"'''''*' 

^ m',l'&[N] 



/ / _ 7'1 «-J2t 



9"^^ P^^ 



['m — m,l — l]e-'^^'^" ^^ 



(m\ 



dl 



where the periodicity of the summand with respect to m' and /' has been exploited in the steps labeled with 
(*). Using (11121 ) in (lllll ) then gives 



triTWrvTWrvl- Ri V If/ ,1 

m,«e[Af] 



(113) 



3j Second term in (I102I I.- Next, we consider the second term on the right-hand side of (I102I I. We have 

tr{rxT,^J ^i* trjrxTg} + , trjrxTm } ^1^ trlr^Tg} = trJTgrx}. (114) 

Using Jo,o = ~^I' we obtain further 

tr{rxT^^J*W»tr{TWrxl} 

= ^/iVtrJTWrxJo^o} 

= ^/iV(TWrx,Jo,o> 

no5i 



Gi) 



^ ^ 2_^ dp,q;m,l Jm,l , Jo,0 



(115) 



4) Approximation of hp^q.- Inserting (11131) and (11151) into (11021) . we obtain the following approximation of 



''P,q- 



h 



p,q 



7 ^ Mp,(j;m,z| + ^|'^p,ij;0,( 



m,ie[iV] 



This can be expressed as 



Pi9 ~ 7 ^ \"'p,q;n,k\ + -'"' 
n,feG[Ar] 



^ Z^ "P,?;",' 



n,fce[Af] 



(116) 



where dp^q-n,k is the 2-D DFT of dp^q-rn,i- We have 



Ojn.n: 



•p,q;n,k 



^^"^ 



, p-J^ikm-nl) 
p,q;m,l e 



m,le[N] 



Got) 1 



^ /4[m',/'] 



m',Pe[Af] 



^ J] ix[m-m',/-/']e-^'i^('=™-"') 



m,l(^[N] 



oJ^" V AM Ai/ 



o 1 



-y= ^ /4[m',/']i?x[n,A;]e' 



-i^(fcm'-™«')eJ2-(S^-i^) 
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^Rx[n,k] Y^ /^[rn',/'] e-^'"[(^-^')'"'-("-^p)''. 



m',l'e[N] 



Eo) 



y/NRx[n,k]^MWv[n-pAn,k-qAk], (117) 

where, as before, An = N/AL and Ak = N/AM. Inserting (11171 ) into (11161 ) finally yields 



hp,q ~ N y^ \Rx[n,k]^MW\j[n—pAn,k — qAk]\ 

n,k&[N] 

which is (|68]). 



2 

n,fce[Af] 
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